#include "poly.h"


////////////////////////////////////////////////////////////////
// Z = POLYISZERO(P) Determine whether the ppololynomial P is zero
//
//   Inputs:
//
//      p       The polynomial to check
//
//   Outputs
//
//       z       True if p is zero, false if p is not
////////////////////////////////////////////////////////////////

int polyIsZero(Vector * p) {

    int i;
    
    // Loop through the vector to look for nonzero coefficients
    for (i = 0; i < p->len; i++) {
        if (p->data[i] != -1) // If a nonzero coefficient is found,
            return 0;         // return false
    }

    // All polynomial coefficients are zero, so return true
    return 1;

}

/////////////////////////////////////////////////////////////////////
// P_TRIMMED = POLYTRIMZEROS(P)  Trim leading zeros from a polynomial
//
//   Inputs:
//
//       p       Polynomial to be trimmed
//
//   Outputs:
//
//       p_trimmed  p with leading zeros removed
/////////////////////////////////////////////////////////////////////

Vector * polyTrimZeros(Vector * p) {
    Vector * p_trimmed;
    int i,trim;

    // If the polynomial is already zero, return the zero polynomial [-1]
    if (polyIsZero(p)) {
        p_trimmed = allocate_vector(1);
        p_trimmed->data[0] = -1;
    }
    else {
        
        //Loop from the end of the vector and count back until a nonzero value
        //is reached.
        trim = 0;
        i = (int)(p->len) - 1;
        while (i >= 0 && p->data[i] == -1) {
            trim++;
            i--;
        }
        
        // Copy the nonzero part of the vector
        p_trimmed = allocate_vector((unsigned long)((int)(p->len) - trim));
        memcpy(p_trimmed->data, p->data, sizeof(int) * ((int)(p->len) - trim));
       
    }
    
    // Return the trimmed polynomial
    return p_trimmed;
    
}

////////////////////////////////////////////////////////
// z = POLYPLUGIN(p.x) Evaluate the polynomial p(x) at x
//
//   Inputs:
//
//       p    The polynomial to be evaluated
//       x    The argument to the polynomial
//
//   Outputs:
//
//       z    The polynomial p(x) evaluated at x
////////////////////////////////////////////////////////
int polyPlugIn(int m, Matrix * add_one, Vector * p, int x) {

    // If the polynomial is zero, return zero
    if (polyIsZero(p)) {
        return -1;
    }
    // If the polynomial only has a constant term, return the constant
    else if (p->len == 0) {
        return p->data[0];
    }
    else {
        
        int i, px, term;
        px = p->data[0];
        
        for (i = 1; i < p->len; i++) {
            if (p->data[i] != -1) {
                term = elementPow(m, add_one, x, i);
                term = elementMult(m, add_one, term, p->data[i]);
                px = elementAdd(m, add_one, px, term);
            }       
        }
        return px;
    }
}


///////////////////////////////////////////////////
// Z = POLYADD(X,Y) Add the two polynomials X and Y
//
//   Inputs:
//
//       x,y     Polynomial summands
//
//   Outputs:
//
//       z       The sum of x and y
///////////////////////////////////////////////////
Vector * polyAdd(int m, Matrix * add_one, Vector * x, Vector * y) {

    Vector *shortPoly, *longPoly, *polySum, *polySumTrimmed;

    // if x is zero, return y
    if (polyIsZero(x)) {
        return copy_vector(y);
    }

    // if y is zero, return x
    else if (polyIsZero(y)) {
        return copy_vector(x);
    }

    else {

        // Determine which polynomial is shorter
        if (x->len > y->len) {
            shortPoly = y;
            longPoly = x;
        }
        else {
            shortPoly = x;
            longPoly = y;
        }

        // Initialize output polynomial vector
        int i;
        polySum = allocate_vector(longPoly->len);
        memcpy(polySum->data, longPoly->data, sizeof(int)*(longPoly->len));

        // Add all coefficients
        for (i = 0; i < shortPoly->len; i++) {
            polySum->data[i] = elementAdd(m, add_one, shortPoly->data[i], longPoly->data[i]);
        }
        
        // Get rid of any leading zeros
        polySumTrimmed = polyTrimZeros(polySum);
        free_vector(polySum);
        
        return polySumTrimmed;
    }
}


/////////////////////////////////////////////////////////
// Z = POLYMULT(X,Y) Multiply the two polynomials X and Y
//
//   Inputs:
//
//       x,y     Polynomial multiplicands
//
//   Outputs:
//
//       z       The product of x and y
/////////////////////////////////////////////////////////
Vector * polyMult(int m, Matrix * add_one, Vector * x, Vector * y) {

    Vector *temp;
 
    // If either polynomial is zero, return zero
    if (polyIsZero(x)) {
        temp = allocate_vector(1);
        temp->data[0] = -1;
        return temp;
    }
    else if (polyIsZero(y)) {
        temp = allocate_vector(1);
        temp->data[0] = -1;
        return temp;
    }

    // Determine which polynomial is shorter
    Vector *shortPoly, *longPoly;
    
    if (x->len > y->len) {
        shortPoly = y;
        longPoly = x;
    }
    else {
        shortPoly = x;
        longPoly = y;
    }

    // Initialize output vector
    Vector * polyProduct;
    
    polyProduct = allocate_vector(x->len + y->len - 1);
    
    // Initialize the vector (can't use memset because we need to initialize 
    // to -1
    int k;
    for (k = 0; k < polyProduct->len; k++) {
        polyProduct->data[k] = -1;
    }

    // Multiply the polynomials.  This is the same operation as convolution.
    int i,j,z;
    Vector * coef;
    Vector * vectZ;

    for (i = 0; i < shortPoly->len; i++) {
        if (i > 0) {
            // Pad by the power of x
            coef = allocate_vector(i);

            for (k = 0; k < coef->len; k++) {
                coef->data[k] = -1;
            }
        }

        // Multiply the coefficients
        if (shortPoly->data[i] != -1) {

            vectZ = allocate_vector(1);

            for (j = 0; j < longPoly->len; j++) {
                
                z = elementMult(m, add_one, shortPoly->data[i], longPoly->data[j]);
                vectZ->data[0] = z;
                if (i > 0 || j > 0) {
                    temp = vector_cat(coef, vectZ);
                    free_vector(coef);
                    coef = temp;
                }
                else {
		    if (i > 0) {
			free_vector(coef);
		    }
                    coef = copy_vector(vectZ);
                }
            }

            free_vector(vectZ);

            temp = polyAdd(m, add_one, polyProduct, coef);
            free_vector(polyProduct);
            polyProduct = temp;

        }   
        if (i > 0 || shortPoly->data[i] != -1) {
            free_vector(coef);
        }
    }
    
    // Get rid of any leading zeros
    temp = polyTrimZeros(polyProduct);
    free_vector(polyProduct);
    polyProduct = temp;

    return polyProduct;
}

/////////////////////////////////////////////////////////////////////
// [Q R] = POLYDIV(N,D) Evaluate the quotient N / D
//
//  Inputs:
//
//       n    Polynomial numerator
//       d    Polynomial denominator
//
//   Outputs:
//
//       q    The quotient n / d
//       r    The remainder of n / d
//
//   Notes:
//
//       Uses pseudocode taken from 
//       http://rosettacode.org/wiki/Polynomial_long_division#Octave
/////////////////////////////////////////////////////////////////////
Vector** polyDiv(int m, Matrix * add_one, Vector *n, Vector *d){
   Vector *q, *r, *dd, *v, *tempVec;
   int i, temp;
   Vector *ntemp;
   Vector **retVect;

   // Check if denominator is zero and print error if so
   if (polyIsZero(d)){
      q = allocate_vector(1);
      q->data[0] = 0;
      r = allocate_vector(1);
      r->data[0] = 0;		// If numerator is 0 return -1 -1
      retVect = (Vector**)malloc(sizeof(Vector*)*2);
      retVect[0] = q;
      retVect[1] = r;
      printf("ERROR: DIVIDE BY ZERO IN POLYDIV!\n");
      return retVect;
   }
   
   // Check if numerator is zero, return zero if true
   if(polyIsZero(n)){
     q = allocate_vector(1);
     q->data[0] = -1;
     r = allocate_vector(1);
     r->data[0] = -1;		// If numerator is 0 return -1 -1
     retVect = (Vector**)malloc(sizeof(Vector*)*2);
     retVect[0] = q;
     retVect[1] = r;
     return retVect;
   }

   // Divide n by d
   ntemp = copy_vector(n); // Use a trash copy of n since n gets changed 
   if (deg(ntemp) >= deg(d)){
      // Create a quotient vector the length of numerator vector
      q = allocate_vector(deg(ntemp)+1);
      // Fill the vector with -1
      for (i=0;i<deg(ntemp)+1; i++){
         q->data[i]=-1;
      }
      // While the degree of the numerator is greater than degree of denominator
      while(deg(ntemp) >= deg(d)){
         // temp holds length of dd vector
         temp = deg(ntemp)-deg(d) + d->len;
         // Allocate a temporary vector dd
         dd = allocate_vector(temp);
         // Fill vector dd
	 for(i=0;i<(deg(ntemp)-deg(d)); i++){
           dd->data[i] = -1;
         }
         for(i=(deg(ntemp)-deg(d)); i< temp;i++){
           dd->data[i] = d->data[i-(deg(ntemp)-deg(d))];
         }
         q->data[deg(ntemp)-deg(d)] = elementDiv(m, add_one, ntemp->data[deg(ntemp)], dd->data[deg(dd)]);
	 v = allocate_vector(1);
         v->data[0] = q->data[deg(ntemp)-deg(d)];
         tempVec = polyMult(m, add_one, dd, v);
         free_vector(dd);
         dd = tempVec;
	 free_vector(v);
         tempVec = polyAdd(m, add_one, ntemp, dd);
         free_vector(ntemp);
         ntemp = tempVec;
	 free_vector(dd);
      }
      r = allocate_vector(ntemp->len);
      free(r->data);
      r->data = ntemp->data;
   }
   else{
      q = allocate_vector(1);      
      q->data[0] = -1;
      r = allocate_vector(ntemp->len);
      free(r->data);
      r->data = ntemp->data;
   }

   tempVec = polyTrimZeros(q);
   free_vector(q);
   q = tempVec;

   tempVec = polyTrimZeros(r);
   free_vector(r);
   r = tempVec;

   retVect = (Vector**)malloc(2*sizeof(Vector*));
   retVect[0] = q;
   retVect[1] = r;

   // Free ntemp.  Note that this is NOT a free_vector() call because ntemp's
   // data was already freed in the call to free r.
   free(ntemp);

   return retVect;
}  

// Helper function to determine the degree of a polynomial
int deg(Vector *p){

  Vector * temp;
  unsigned long len;

  if (polyIsZero(p)){
    return -1;
  }
  temp = polyTrimZeros(p);
  len = temp->len;
  free_vector(temp);
  
  return (len - 1);
  
}

// findRoots function
Vector* findRoots(int m, Matrix *add_one, Vector *p){
  Vector *rts;
  int i, px;
  unsigned long temp;
  int *newdata;
  int root_counter;
  int rtsFull = 0;	// Variable to determine if rts has been created yet
  
  root_counter = 0;

  // Loop until vector is made
  for(i=0;i<(1<<m); i++){
     // Determine if result is 0
     px = polyPlugIn(m, add_one, p,i-1);

     if (px == -1) root_counter++;
     // If plug in result is 0 add it to root (rts) vector
     if ((px == -1) && (rtsFull == 0)){	// Case where rts not yet made
        rts = allocate_vector(1);
        rts->data[0] = i-1;
        rtsFull = 1;
     }
     else if (px == -1){ // Case where rts vector is expanding
        temp = rts->len + 1;
    	newdata = malloc((rts->len)*sizeof(int));
        memcpy(newdata, rts->data, (rts->len) * sizeof(int));
        free_vector(rts);
        rts = allocate_vector(temp);
        for(temp = 0; temp<rts->len-1; temp++){
           rts->data[temp] = newdata[temp];
        }
        rts->data[rts->len-1] = i-1;
        free(newdata);
     }       
  }

  if (root_counter > 0) {
      return rts; // Return root (rts) vector 
  }
  else {
      return NULL;
  }
}
